Analisis Descriptivo de Series de Tiempo en R

Gráfico de una Serie de tiempo

Vamos a analizar de forma descriptiva algunas serie de tiempo. Empezaremos por la serie de desempleo de los Estados Unidos que viene en el paquete TSstudio.

library(TSstudio)
data(USUnRate)
ts_info(USUnRate)

##  The USUnRate series is a ts object with 1 variable and 864 observations
##  Frequency: 12 
##  Start time: 1948 1 
##  End time: 2019 12

class(USUnRate)

## [1] "ts"

plot(USUnRate,main = "US Monthly Unemployment Rate",ylab="Unemployment Rate (%)")

La serie es medida mensual, es decir, presenta una frecuencia de 12. Qué características podemos observar? * Tendencia? * Estacionalidad? * Cíclos? * Varianza marginal no constante? Vamos a seleccionar un periodo de tiempo mas corto

unemployment <- window(USUnRate, start = c(1990,1))
ts_plot(unemployment,
          title = "US Monthly Unemployment Rate",
           Ytitle = "Unemployment Rate (%)",
           Xtitle = "Year",
          Xgrid = TRUE,
Ygrid = TRUE)

Note que aquí podemos ver varias características: Estacionalidad(no es tan evidente), tres periodos de ciclos, el primero de 1990 a 2000,el segundo de 2000 a 2007, y el tercero de 2007 a 2019. No parece tener una heterocedasticidad marginal.

Veamos ahora la tasa de desempleo de Colombia. Hay que hacerle un ajuste a la base de datos porque está en orden descendente en el tiempo.

Desempleo y Empleo Colombia

library(readxl)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

DesempleoyEmpleo <- read_excel("DesempleoyEmpleo.xlsx", range="A6:C269")
str(DesempleoyEmpleo)

## tibble [263 × 3] (S3: tbl_df/tbl/data.frame)
##  $ Fecha            : chr [1:263] "2022-11" "2022-10" "2022-09" "2022-08" ...
##  $ Tasa_de_ocupacion: num [1:263] 57.4 57.7 57.2 56.7 56.5 ...
##  $ Tasa_de_desempleo: num [1:263] 9.5 9.72 10.75 10.63 10.99 ...

DesempleoyEmpleo_1=DesempleoyEmpleo %>% map_df(rev)
tail(DesempleoyEmpleo)

## # A tibble: 6 × 3
##   Fecha   Tasa_de_ocupacion Tasa_de_desempleo
##   <chr>               <dbl>             <dbl>
## 1 2001-06              56.2              15.3
## 2 2001-05              56.2              14.0
## 3 2001-04              55.8              14.5
## 4 2001-03              57.6              15.8
## 5 2001-02              56.9              17.4
## 6 2001-01              57.6              16.6

head(DesempleoyEmpleo_1)

## # A tibble: 6 × 3
##   Fecha   Tasa_de_ocupacion Tasa_de_desempleo
##   <chr>               <dbl>             <dbl>
## 1 2001-01              57.6              16.6
## 2 2001-02              56.9              17.4
## 3 2001-03              57.6              15.8
## 4 2001-04              55.8              14.5
## 5 2001-05              56.2              14.0
## 6 2001-06              56.2              15.3

library(zoo)

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(xts)

## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################

## 
## Attaching package: 'xts'

## The following objects are masked from 'package:dplyr':
## 
##     first, last

Fechas=as.yearmon(DesempleoyEmpleo_1$Fecha)
Desempleo_Col_xts=xts(x = DesempleoyEmpleo_1$Tasa_de_desempleo,frequency = 12,order.by = Fechas)
ts_info(Desempleo_Col_xts)

##  The Desempleo_Col_xts series is a xts object with 1 variable and 263 observations
##  Frequency: monthly 
##  Start time: Jan 2001 
##  End time: Nov 2022

plot(Desempleo_Col_xts)

ts_plot(Desempleo_Col_xts,
           title = "Tasa de Desemplo Mensual Colombia",
           Ytitle = "Tasa de Desempleo(%)",
           Xtitle = "Año",
           Xgrid = TRUE,
Ygrid = TRUE)

Qué características presenta esta serie?

Otros objetos

Note que también se puede crear un objeto tsibble y usar el paquete feast, fable y fabletools(tsibble) y timetk(tibble), lo cuales funcionan con el paquete tidyverse.

require(feasts)

## Loading required package: feasts

## Loading required package: fabletools

require(fable)

## Loading required package: fable

require(timetk)

## Loading required package: timetk

require(tsibble)

## Loading required package: tsibble

## 
## Attaching package: 'tsibble'

## The following object is masked from 'package:zoo':
## 
##     index

## The following object is masked from 'package:lubridate':
## 
##     interval

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union

require(lubridate)
###Creación objeto tssible a partir de un objeto tibble
df_desempleo=data.frame(Desempleo=DesempleoyEmpleo_1$Tasa_de_desempleo,Fecha=DesempleoyEmpleo_1$Fecha)
tbl_desempleo=tibble(df_desempleo)
tbl_desempleo_format_fecha=tbl_desempleo
tbl_desempleo_format_fecha$Fecha=yearmonth(tbl_desempleo_format_fecha$Fecha)
###El tipo de fechas debe ser alguno que reconozca tsibble

tsbl_desempleo=as_tsibble(tbl_desempleo_format_fecha,index=Fecha)   ####La fecha en tsibble es importante

##Gráfica de tsibble
autoplot(tsbl_desempleo,Desempleo)+labs(tittle="Serie de Desempleo Colombia Mensual",y="Tasa de Desempleo")

tbl_desempleo$Fecha<-as.Date(zoo::as.yearmon(tbl_desempleo$Fecha))
tbl_desempleo

## # A tibble: 263 × 2
##    Desempleo Fecha     
##        <dbl> <date>    
##  1      16.6 2001-01-01
##  2      17.4 2001-02-01
##  3      15.8 2001-03-01
##  4      14.5 2001-04-01
##  5      14.0 2001-05-01
##  6      15.3 2001-06-01
##  7      15.2 2001-07-01
##  8      14.9 2001-08-01
##  9      14.1 2001-09-01
## 10      14.5 2001-10-01
## # ℹ 253 more rows

###Gráfica timetk
tbl_desempleo%>%plot_time_series(.value=Desempleo,.date_var=Fecha)

Análisis de Tendencias

Vamos a ver la forma de estimar la tendencia y/o eliminarla.

library(astsa)
library(TSstudio)
data(chicken)
ts_info(chicken)

##  The chicken series is a ts object with 1 variable and 180 observations
##  Frequency: 12 
##  Start time: 2001 8 
##  End time: 2016 7

plot(chicken,main="Precio Mensual de la Libra de Pollo en Estados Unidos", ylab="Precio en Centavos de Dólar")

#ts_plot(chicken)

Al parecer la serie de precios mensuales del pollo presenta una tendencia creciente al parecer lineal, es decir \[y_t=\mu_t+a_t\] o mas específicamente

\[y_t=\beta_0+\beta_1 t +a_t\]

Un primer intento de ver como sería la tendencia sería a través del uso del filtro de promedios móviles.

Descomposición por Filtro de Promedios Móviles

Vamos a usar la descomposición por medio de filtros de promedios móviles

chicken_decompo=decompose(chicken)
plot(chicken_decompo)

chicken_decompo$trend

##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2001                                                                      
## 2002        NA  63.91958  63.75667  63.54250  63.32875  63.15500  63.05458
## 2003  63.64792  63.96792  64.36708  64.81542  65.31667  65.89792  66.51458
## 2004  71.92417  72.96750  73.82292  74.50458  75.07875  75.53583  75.88917
## 2005  75.35458  74.87458  74.52875  74.33750  74.18583  74.00208  73.75333
## 2006  71.02208  70.63875  70.27000  69.88542  69.53458  69.30333  69.28708
## 2007  73.70708  74.62875  75.53333  76.40667  77.19292  77.87083  78.43000
## 2008  80.89500  81.48792  82.07125  82.68125  83.38750  84.19292  85.03333
## 2009  87.24208  87.18625  86.97083  86.62875  86.23333  85.83042  85.45208
## 2010  84.68500  84.69750  84.85958  85.14083  85.44125  85.71333  85.92833
## 2011  86.39167  86.38500  86.45042  86.59625  86.84792  87.19125  87.60042
## 2012  91.02625  91.62125  92.18625  92.74917  93.34292  93.97792  94.66958
## 2013  99.52667 100.49167 101.40917 102.23292 102.95292 103.56333 104.05833
## 2014 106.44125 106.96125 107.52875 108.20167 108.95417 109.73583 110.53667
## 2015 114.26083 114.50958 114.68250 114.76000 114.75208 114.69875 114.60333
## 2016 113.03167        NA        NA        NA        NA        NA        NA
##            Aug       Sep       Oct       Nov       Dec
## 2001        NA        NA        NA        NA        NA
## 2002  63.03542  63.09125  63.18125  63.27625  63.41958
## 2003  67.17167  67.90875  68.76083  69.72792  70.79583
## 2004  76.14000  76.26292  76.26458  76.13750  75.82708
## 2005  73.41375  72.99042  72.48750  71.95000  71.45333
## 2006  69.53958  70.06750  70.84500  71.77125  72.74708
## 2007  78.90083  79.32750  79.69417  80.02250  80.39333
## 2008  85.76458  86.26667  86.59333  86.87833  87.12667
## 2009  85.13500  84.92125  84.84500  84.81958  84.75667
## 2010  86.08375  86.24417  86.37750  86.42792  86.42208
## 2011  88.07750  88.61125  89.17625  89.77500  90.40333
## 2012  95.41000  96.14625  96.89542  97.70167  98.58000
## 2013 104.45875 104.79708 105.15125 105.54292 105.96083
## 2014 111.32708 112.08917 112.78208 113.39792 113.91000
## 2015 114.46792 114.28542 114.03375 113.72917 113.39000
## 2016

Simule unos datos III, crear un objeto con periodicidad s=12, y extraer la descomposición, qué puede observar?

y=ts(rnorm(1000,0,1),start=c(2000,01),frequency = 12)
plot(decompose(y))

Se puede calcular manualmente los filtros

chicken_decompo$trend

##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2001                                                                      
## 2002        NA  63.91958  63.75667  63.54250  63.32875  63.15500  63.05458
## 2003  63.64792  63.96792  64.36708  64.81542  65.31667  65.89792  66.51458
## 2004  71.92417  72.96750  73.82292  74.50458  75.07875  75.53583  75.88917
## 2005  75.35458  74.87458  74.52875  74.33750  74.18583  74.00208  73.75333
## 2006  71.02208  70.63875  70.27000  69.88542  69.53458  69.30333  69.28708
## 2007  73.70708  74.62875  75.53333  76.40667  77.19292  77.87083  78.43000
## 2008  80.89500  81.48792  82.07125  82.68125  83.38750  84.19292  85.03333
## 2009  87.24208  87.18625  86.97083  86.62875  86.23333  85.83042  85.45208
## 2010  84.68500  84.69750  84.85958  85.14083  85.44125  85.71333  85.92833
## 2011  86.39167  86.38500  86.45042  86.59625  86.84792  87.19125  87.60042
## 2012  91.02625  91.62125  92.18625  92.74917  93.34292  93.97792  94.66958
## 2013  99.52667 100.49167 101.40917 102.23292 102.95292 103.56333 104.05833
## 2014 106.44125 106.96125 107.52875 108.20167 108.95417 109.73583 110.53667
## 2015 114.26083 114.50958 114.68250 114.76000 114.75208 114.69875 114.60333
## 2016 113.03167        NA        NA        NA        NA        NA        NA
##            Aug       Sep       Oct       Nov       Dec
## 2001        NA        NA        NA        NA        NA
## 2002  63.03542  63.09125  63.18125  63.27625  63.41958
## 2003  67.17167  67.90875  68.76083  69.72792  70.79583
## 2004  76.14000  76.26292  76.26458  76.13750  75.82708
## 2005  73.41375  72.99042  72.48750  71.95000  71.45333
## 2006  69.53958  70.06750  70.84500  71.77125  72.74708
## 2007  78.90083  79.32750  79.69417  80.02250  80.39333
## 2008  85.76458  86.26667  86.59333  86.87833  87.12667
## 2009  85.13500  84.92125  84.84500  84.81958  84.75667
## 2010  86.08375  86.24417  86.37750  86.42792  86.42208
## 2011  88.07750  88.61125  89.17625  89.77500  90.40333
## 2012  95.41000  96.14625  96.89542  97.70167  98.58000
## 2013 104.45875 104.79708 105.15125 105.54292 105.96083
## 2014 111.32708 112.08917 112.78208 113.39792 113.91000
## 2015 114.46792 114.28542 114.03375 113.72917 113.39000
## 2016

as.numeric(chicken_decompo$trend)

##   [1]        NA        NA        NA        NA        NA        NA  63.91958
##   [8]  63.75667  63.54250  63.32875  63.15500  63.05458  63.03542  63.09125
##  [15]  63.18125  63.27625  63.41958  63.64792  63.96792  64.36708  64.81542
##  [22]  65.31667  65.89792  66.51458  67.17167  67.90875  68.76083  69.72792
##  [29]  70.79583  71.92417  72.96750  73.82292  74.50458  75.07875  75.53583
##  [36]  75.88917  76.14000  76.26292  76.26458  76.13750  75.82708  75.35458
##  [43]  74.87458  74.52875  74.33750  74.18583  74.00208  73.75333  73.41375
##  [50]  72.99042  72.48750  71.95000  71.45333  71.02208  70.63875  70.27000
##  [57]  69.88542  69.53458  69.30333  69.28708  69.53958  70.06750  70.84500
##  [64]  71.77125  72.74708  73.70708  74.62875  75.53333  76.40667  77.19292
##  [71]  77.87083  78.43000  78.90083  79.32750  79.69417  80.02250  80.39333
##  [78]  80.89500  81.48792  82.07125  82.68125  83.38750  84.19292  85.03333
##  [85]  85.76458  86.26667  86.59333  86.87833  87.12667  87.24208  87.18625
##  [92]  86.97083  86.62875  86.23333  85.83042  85.45208  85.13500  84.92125
##  [99]  84.84500  84.81958  84.75667  84.68500  84.69750  84.85958  85.14083
## [106]  85.44125  85.71333  85.92833  86.08375  86.24417  86.37750  86.42792
## [113]  86.42208  86.39167  86.38500  86.45042  86.59625  86.84792  87.19125
## [120]  87.60042  88.07750  88.61125  89.17625  89.77500  90.40333  91.02625
## [127]  91.62125  92.18625  92.74917  93.34292  93.97792  94.66958  95.41000
## [134]  96.14625  96.89542  97.70167  98.58000  99.52667 100.49167 101.40917
## [141] 102.23292 102.95292 103.56333 104.05833 104.45875 104.79708 105.15125
## [148] 105.54292 105.96083 106.44125 106.96125 107.52875 108.20167 108.95417
## [155] 109.73583 110.53667 111.32708 112.08917 112.78208 113.39792 113.91000
## [162] 114.26083 114.50958 114.68250 114.76000 114.75208 114.69875 114.60333
## [169] 114.46792 114.28542 114.03375 113.72917 113.39000 113.03167        NA
## [176]        NA        NA        NA        NA        NA

m_filter=12
filter_1=stats::filter(chicken, filter = rep(1/m_filter, m_filter), sides = 2)
filter_2=forecast::ma(chicken,m_filter)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

filter_3=forecast::ma(chicken,m_filter,centre = T)

Tarea: Estimar yEliminar la tendencia por un filtro de promedio móviles, y a la serie resultante extraerle la estacionalidad a partir de otro filtro de promedios móviles.(note que para este caso lo que hay que hacer es agrupar por meses a través de todos los años y calcular el promedio.)

Para mas detalles https://otexts.com/fpp3/moving-averages.html.

Ajuste de tedencia determinística

summary(fit <- lm(chicken~time(chicken), na.action=NULL))

## 
## Call:
## lm(formula = chicken ~ time(chicken), na.action = NULL)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.7411 -3.4730  0.8251  2.7738 11.5804 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -7.131e+03  1.624e+02  -43.91   <2e-16 ***
## time(chicken)  3.592e+00  8.084e-02   44.43   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.696 on 178 degrees of freedom
## Multiple R-squared:  0.9173, Adjusted R-squared:  0.9168 
## F-statistic:  1974 on 1 and 178 DF,  p-value: < 2.2e-16

plot(chicken, ylab="centavos por libra") 
abline(fit,col = "red") # Se añade la recta ajusta

###Eliminamos la tendencia con la predicción la recta
ElimiTendchick=chicken-predict(fit)
plot(ElimiTendchick,main="Serie Chicken Sin tendencia")

acf(ElimiTendchick,lag.max =179 )

library(tidyverse)
library(lubridate)
library(timetk)
library(tsibble)

# Se configura para gráficas plotly(# FALSE retorna ggplots y no plotly)
interactive <- FALSE
indice_chicken=as.Date(as.yearmon(tk_index(chicken)))
## Otra forma de extraer el indice estimetk::tk_index(chicken)
df_chicken=data.frame(Fecha=indice_chicken,Pollo=as.matrix(chicken))
str(df_chicken)

## 'data.frame':    180 obs. of  2 variables:
##  $ Fecha: Date, format: "2001-08-01" "2001-09-01" ...
##  $ Pollo: num  65.6 66.5 65.7 64.3 63.2 ...

tibble_chicken=tibble(df_chicken)
duplicates(tibble_chicken, key = NULL, index=Fecha)   ##Mirar si hay registros duplicados

## # A tibble: 0 × 2
## # ℹ 2 variables: Fecha <date>, Pollo <dbl>

tibble_chicken_fechas_correct=tibble_chicken
tibble_chicken_fechas_correct$Fecha=yearmonth(tibble_chicken_fechas_correct$Fecha)
print(duplicates(tibble_chicken_fechas_correct, key = NULL, index=Fecha))

## # A tibble: 0 × 2
## # ℹ 2 variables: Fecha <mth>, Pollo <dbl>

tsibble_chicken=tsibble(tibble_chicken_fechas_correct,index=Fecha)

tsibble_chicken

## # A tsibble: 180 x 2 [1M]
##       Fecha Pollo
##       <mth> <dbl>
##  1 2001 Aug  65.6
##  2 2001 Sep  66.5
##  3 2001 Oct  65.7
##  4 2001 Nov  64.3
##  5 2001 Dec  63.2
##  6 2002 Jan  62.9
##  7 2002 Feb  62.9
##  8 2002 Mar  62.7
##  9 2002 Apr  62.5
## 10 2002 May  63.4
## # ℹ 170 more rows

Note que se ajusta una tendencia que se ve que no es lineal.

tibble_chicken%>%timetk::plot_time_series(Fecha, Pollo, 
                   .interactive = TRUE,
                   .plotly_slider = TRUE)

###Usa Loess para hacer el ajuste de la tendencia, es decir usar smooth_vec() como versión simplificada de stats::loess()
tibble_chicken%>%mutate(Pollo_ajus=smooth_vec(Pollo,span = 0.35, degree = 2))

## # A tibble: 180 × 3
##    Fecha      Pollo Pollo_ajus
##    <date>     <dbl>      <dbl>
##  1 2001-08-01  65.6       62.8
##  2 2001-09-01  66.5       62.8
##  3 2001-10-01  65.7       62.9
##  4 2001-11-01  64.3       63.0
##  5 2001-12-01  63.2       63.2
##  6 2002-01-01  62.9       63.3
##  7 2002-02-01  62.9       63.5
##  8 2002-03-01  62.7       63.6
##  9 2002-04-01  62.5       63.8
## 10 2002-05-01  63.4       64.0
## # ℹ 170 more rows

tibble_chicken%>%mutate(Pollo_ajus=smooth_vec(Pollo,span = 0.75, degree = 2))%>%
  ggplot(aes(Fecha, Pollo)) +
    geom_line() +
    geom_line(aes(y = Pollo_ajus), color = "red")

Descomposición STL

STL son las iniciales de “Seasonal and Trend decomposition using Loess”,el cual fue desarrollado por R. B. Cleveland et al. (1990).

library(feasts)
library(fable)

tsibble_chicken<-as_tsibble(chicken)
str(tsibble_chicken)

## tbl_ts [180 × 2] (S3: tbl_ts/tbl_df/tbl/data.frame)
##  $ index: mth [1:180] 2001 Aug, 2001 Sep, 2001 Oct, 2001 Nov, 2001 Dec, 2002 Jan...
##  $ value: num [1:180] 65.6 66.5 65.7 64.3 63.2 ...
##  - attr(*, "key")= tibble [1 × 1] (S3: tbl_df/tbl/data.frame)
##   ..$ .rows: list<int> [1:1] 
##   .. ..$ : int [1:180] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..@ ptype: int(0) 
##  - attr(*, "index")= chr "index"
##   ..- attr(*, "ordered")= logi TRUE
##  - attr(*, "index2")= chr "index"
##  - attr(*, "interval")= interval [1:1] 1M
##   ..@ .regular: logi TRUE

tsibble_chicken %>%
  model(
    STL(value ~ trend() +
                   season(window = "periodic"),
    robust = TRUE)) %>%
  components() %>%
  autoplot()

Note que en ambos casos se obliga a extraer un componente estacional, sin embargo puede que está componente en verdad no exista, por eso se debe verificar que en efecto hay.

set.seed(154) 
w = rnorm(200); x = cumsum(w) 
wd = w +.2; xd = cumsum(wd)
plot.ts(xd, ylim=c(-5,55), main="Caminata Aletoria", ylab='')
lines(x, col=4); abline(h=0, col=4, lty=2); abline(a=0, b=.2, lty=2)

La diferencia ordinaria de orden 1 es , \[\nabla^1 Y_t=(1-B)^1 Y_t=Y_t-Y_{t-1}\]

dx=diff(x)
plot.ts(dx, ylim=c(-5,5), main="Serie Diferenciada", ylab='')

####Con drift
dxd=diff(xd)
plot.ts(dxd, ylim=c(-5,5), main="Serie con drift Diferenciada", ylab='')

par(mar = c(2,2,2,2))
fit = lm(chicken~time(chicken), na.action=NULL) # Regresión sobre el tiempo
par(mfrow=c(2,1))
plot(resid(fit), type="l", main="sin tendencia") 
plot(diff(chicken), type="l", main="Primera Diferencia") 

par(mar = c(3,2,3,2))
par(mfrow=c(3,1)) # plot ACFs
acf(chicken, 48, main="ACF Pollo")
acf(resid(fit), 48, main="ACF Sin tendencia") 
acf(diff(chicken), 48, main="ACF Primera Diferencia")

Transformación de Box-Cox para Estabilizar la Varianza Marginal

En ocasiones la serie presenta varianza marginal no constante a lo largo del tiempo, lo cual hace necesario tener en cuenta tal característica. En este caso, se siguiere hacer una transformación de potencia para estabilizar la varianza. Esta familia de transformaciones se llaman transformaciones Box-Cox.

\[ f_{\lambda}(u_{t})= \begin{cases} \lambda^{-1}(u^{\lambda}_{t}-1), & \text{si $u_{t} \geq 0$, para $\lambda>0$,}\\ \ln(u_{t}), &\text{ si $u_{t}>0$, para $\lambda=0$}. \end{cases} \]

data("AirPassengers")
plot(AirPassengers)

#####Transformación Box-Cox
#library(FitAR)
library(forecast)

## 
## Attaching package: 'forecast'

## The following object is masked from 'package:astsa':
## 
##     gas

## The following object is masked from 'package:fabletools':
## 
##     accuracy

forecast::BoxCox.lambda(AirPassengers, method ="loglik", lower = -1, upper = 3) ###Me entrega el valor de lambda 

## [1] 0.2

##method="loglik"
#FitAR::BoxCox(AirPassengers)###Me entrega una gráfica
plot(forecast::BoxCox(AirPassengers,lambda=0.2))###

lAirPass=log(AirPassengers)
#x11()
par(mar = c(1,1,1,1))
par(mfrow=c(2,1))
plot(AirPassengers,main="Serie de Pasajeros sin Transformar")
plot(lAirPass,main="Series con Transformación BoxCox")

forecast::BoxCox.lambda(lAirPass, method ="guerrero", lower = -1, upper = 3)

## [1] -0.4190648

##Box-Cox con timetk
timetk::box_cox_vec(AirPassengers,lambda = 'auto',silent = F)

## box_cox_vec(): Using value for lambda: -0.294715585559316

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1949 2.548484 2.561375 2.588408 2.582937 2.567506 2.593720 2.615089 2.615089
## 1950 2.555038 2.577299 2.603899 2.593720 2.575381 2.616631 2.646225 2.646225
## 1951 2.610379 2.618160 2.656279 2.636912 2.648795 2.656279 2.680103 2.680103
## 1952 2.647515 2.658701 2.673640 2.659900 2.662270 2.699009 2.709884 2.720049
## 1953 2.676904 2.676904 2.715050 2.714201 2.709007 2.720866 2.737089 2.742835
## 1954 2.685298 2.668053 2.714201 2.707236 2.713347 2.737089 2.762580 2.756933
## 1955 2.720049 2.712489 2.739270 2.740706 2.741419 2.770363 2.796341 2.787869
## 1956 2.751056 2.746317 2.771524 2.769193 2.772100 2.801088 2.818144 2.814820
## 1957 2.770363 2.761963 2.792420 2.788382 2.791921 2.821786 2.837892 2.838594
## 1958 2.784223 2.772100 2.795371 2.788382 2.795857 2.826872 2.846724 2.851232
## 1959 2.794394 2.785275 2.815241 2.810978 2.820985 2.840332 2.864126 2.867216
## 1960 2.819775 2.808794 2.820583 2.836477 2.840332 2.860370 2.883509 2.879580
##           Sep      Oct      Nov      Dec
## 1949 2.595457 2.563441 2.529834 2.561375
## 1950 2.629937 2.590196 2.552878 2.602242
## 1951 2.663443 2.635540 2.611963 2.640966
## 1952 2.690331 2.671428 2.648795 2.674735
## 1953 2.715895 2.692301 2.658701 2.682201
## 1954 2.733382 2.709007 2.684272 2.709007
## 1955 2.768604 2.744238 2.715895 2.747003
## 1956 2.791921 2.765020 2.742129 2.765020
## 1957 2.814399 2.787869 2.764414 2.782096
## 1958 2.814399 2.793903 2.767420 2.782631
## 1959 2.837187 2.815659 2.795371 2.814820
## 1960 2.852177 2.836477 2.808352 2.825716

Note que ahora usamos la misma función para verificar si en verdad la varianza fue estabilizada.

#FitAR::BoxCox(lAirPass)
forecast::BoxCox.lambda(lAirPass, method = "guerrero", lower = -1, upper = 2)

## [1] -0.419081

library(VGAM)

## Loading required package: stats4

## Loading required package: splines

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:VGAM':
## 
##     logit

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

VGAM::yeo.johnson(AirPassengers, lambda = 0)

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1949 4.727388 4.779123 4.890349 4.867534 4.804021 4.912655 5.003946 5.003946
## 1950 4.753590 4.844187 4.955827 4.912655 4.836282 5.010635 5.141664 5.141664
## 1951 4.983607 5.017280 5.187386 5.099866 5.153292 5.187386 5.298317 5.298317
## 1952 5.147494 5.198497 5.267858 5.204007 5.214936 5.389072 5.442418 5.493061
## 1953 5.283204 5.283204 5.468060 5.463832 5.438079 5.497168 5.579730 5.609472
## 1954 5.323010 5.241747 5.463832 5.429346 5.459586 5.579730 5.713733 5.683580
## 1955 5.493061 5.455321 5.590987 5.598422 5.602119 5.755742 5.899897 5.852202
## 1956 5.652489 5.627621 5.762051 5.749393 5.765191 5.926926 6.025866 6.006353
## 1957 5.755742 5.710427 5.877736 5.855072 5.874931 6.047372 6.144186 6.148468
## 1958 5.831882 5.765191 5.894403 5.855072 5.897154 6.077642 6.198479 6.226537
## 1959 5.888878 5.837730 6.008813 5.983936 6.042633 6.159095 6.308098 6.327937
## 1960 6.035481 5.971262 6.040255 6.135565 6.159095 6.284134 6.434547 6.408529
##           Sep      Oct      Nov      Dec
## 1949 4.919981 4.787492 4.653960 4.779123
## 1950 5.068904 4.897840 4.744932 4.948760
## 1951 5.220356 5.093750 4.990433 5.117994
## 1952 5.347108 5.257495 5.153292 5.273000
## 1953 5.472271 5.356586 5.198497 5.308268
## 1954 5.560682 5.438079 5.318120 5.438079
## 1955 5.746203 5.616771 5.472271 5.631212
## 1956 5.874931 5.726848 5.605802 5.726848
## 1957 6.003887 5.852202 5.723585 5.820083
## 1958 6.003887 5.886104 5.739793 5.823046
## 1959 6.139885 6.011267 5.894403 6.006353
## 1960 6.232448 6.135565 5.968708 6.070738

car::yjPower(AirPassengers,lambda=0)

##   [1] 4.727388 4.779123 4.890349 4.867534 4.804021 4.912655 5.003946 5.003946
##   [9] 4.919981 4.787492 4.653960 4.779123 4.753590 4.844187 4.955827 4.912655
##  [17] 4.836282 5.010635 5.141664 5.141664 5.068904 4.897840 4.744932 4.948760
##  [25] 4.983607 5.017280 5.187386 5.099866 5.153292 5.187386 5.298317 5.298317
##  [33] 5.220356 5.093750 4.990433 5.117994 5.147494 5.198497 5.267858 5.204007
##  [41] 5.214936 5.389072 5.442418 5.493061 5.347108 5.257495 5.153292 5.273000
##  [49] 5.283204 5.283204 5.468060 5.463832 5.438079 5.497168 5.579730 5.609472
##  [57] 5.472271 5.356586 5.198497 5.308268 5.323010 5.241747 5.463832 5.429346
##  [65] 5.459586 5.579730 5.713733 5.683580 5.560682 5.438079 5.318120 5.438079
##  [73] 5.493061 5.455321 5.590987 5.598422 5.602119 5.755742 5.899897 5.852202
##  [81] 5.746203 5.616771 5.472271 5.631212 5.652489 5.627621 5.762051 5.749393
##  [89] 5.765191 5.926926 6.025866 6.006353 5.874931 5.726848 5.605802 5.726848
##  [97] 5.755742 5.710427 5.877736 5.855072 5.874931 6.047372 6.144186 6.148468
## [105] 6.003887 5.852202 5.723585 5.820083 5.831882 5.765191 5.894403 5.855072
## [113] 5.897154 6.077642 6.198479 6.226537 6.003887 5.886104 5.739793 5.823046
## [121] 5.888878 5.837730 6.008813 5.983936 6.042633 6.159095 6.308098 6.327937
## [129] 6.139885 6.011267 5.894403 6.006353 6.035481 5.971262 6.040255 6.135565
## [137] 6.159095 6.284134 6.434547 6.408529 6.232448 6.135565 5.968708 6.070738

Ejercicio: Use la serie varve del paquete astsa para chequear si es necesario hacer transformación Box-Cox.

Vamos a correr lo mismo pero en Python para la transformación de BoxCox

#```{r configuración, include=FALSE} #library(knitr) #library(reticulate) ###Nos ubicamos en el ambiente dentro de la terminal y damos which python y copiamos la ruta “/opt/anaconda3/envs/Python38andR/bin/python” ##“/Users/macbook/opt/anaconda3/envs/Ambiente38/bin/python” #/Users/sergiocalderonunal/opt/anaconda3/envs/Ambiente3104/bin/python #use_python(“/Users/macbook/opt/anaconda3/envs/Ambiente38/bin/python”) #use_python(“/Users/sergiocalderonunal/opt/anaconda3/envs/Ambiente3104/bin/python”) #use_virtualenv(“~/Python38andR”) #py_config() #pd <- import(“pandas”) #arreglo<-pd\(array(c(1, 2, 3)) #print(arreglo) #arreglo\)shape

#```

#```{python inicio} #import sys #print(sys.path) #import pandas as pd #import matplotlib.pylab as plt #data = pd.read_csv(‘AirPassengers.csv’) #print(data) #print(‘Data Types:’) #print(data.dtypes) #con=data[‘Month’] #data[‘Month’]=pd.to_datetime(data[‘Month’]) #data #pasajeros=data.set_index(‘Month’) #check datatype of index

#convert to time series: #ts = pasajeros[‘NPassengers’] #ts.head(10)

####Graficar la Serie##### #plt.plot(ts) #plt.title(‘AirPassengers’) #plt.show()

#```

#{python inicio_1} #print(r.AirPassengers) #print(r.lAirPass) #

#```{r uso de objeto creado en python}

#py\(data #library(ggplot2) #ggplot2::ggplot(data = py\)data,aes(x=Month,y=NPassengers) #)+geom_line()

#```

Gráficas de Retardos

Vamos a hacer gráficos de dispersión para chequear que tipos de relaciones hay entre los retardos de la variable interés. Vamos a trabajar con algunas series, por ejemplo: * Indice ambiental mensual (soi)(Southern Oscillation Index), el cual mide los cambios en la presión del aire, relacionados con las temperaturas de la superficie del mar en el Océano Pacífico central. * la serie rec (reclutamiento asociada al soi), número de nuevos peces. * Consumo mensual de gas natural en EE. UU.(USgas) medido en Billones de pies cúbicos. Esto permite chequear si hay posibles relaciones no-lineales.

library(astsa)
data("soi")
ts_info(soi)

##  The soi series is a ts object with 1 variable and 453 observations
##  Frequency: 12 
##  Start time: 1950 1 
##  End time: 1987 9

?soi

par(mar = c(2,2,2,2))
plot(soi, main="Indice soi")

par(mar = c(3,2,3,2))
astsa::lag1.plot(soi, 12,corr=F)  ###El 12 indica cuantos retardos y_t-k contra y_t 

###Hacer la gráfica con x11()

En el gráfico de dispersión podemos ver que se muestra un ajuste no paramétrico, de la posible relación entre las variables al igual que una estimación de la autocorrelación entre \(s_t\) y \(s_{t-h}\). Vemos que varias de las relaciones exploradas parecen ser lineales. Uno puede observar que existen relaciones lineales positivas en los rezagos \(h = 1, 2, 11, 12\), mientras que negativas en los rezagos \(h=6,7\), las demás parecen ser no significativas o no lineales.

#pdf('/Users/sergiocalderon/Documents/Documentos - iMac de Sergio/Documentos iMac Sergio/Notas de Clase/Notas de clase/Notas de Clase Series de Tiempo Univariadas/Graficas/DispersionSoiRec.pdf',paper="USr")
par(mar = c(3,2,3,2))
lag2.plot(soi, rec, 8)   #El 2 de lag2.plot es porque intervenienen dos serie de tiempo.

#dev.off()
lag2.plot(soi, rec, 8,corr=F)

Note también que el gráfico de dispersión del índice de nuevos peces con retardos del indice soi nos muestra también relaciones posiblemente lineales, al igual que no lineales. Ahora, podemos crear el gráfico de autocorrelación simple que nos permite estimar y graficar las autocorrelaciones para diferentes rezagos.

Note que con la función ts_lags de TSstudio podemos hacer una gráfica similar.

ts_lags(soi,lags=1:12)

La función de autocorrelación simple o gráfico ACF y PACF

Cuando el proceso es estacionario, o al menos no presenta tendencia, podemos usar el gráfico acf para explorar las posibles relaciones lineales a diferentes rezagos. En seguida mostramos la función de autocorrelación para el índice soi y la serie de nuevos peces.

par(mfrow=c(2,1))
par(mar = c(2.7,2,2.7,2))
acf(soi, 48, main="Southern Oscillation Index")
acf(rec, 48, main="Recruitment")

#dev.off()

ccf(rec, soi, 48, main="SOI vs Recruitment", ylab="CCF")

#Corr(rec_{t},soi_{t-h})
#Corr(soi_{t},rec_{t-h},)

index_soi_rec=as.Date(as.yearmon(tk_index(soi)))

df_soi_rec=data.frame(Fecha_soi_rec=index_soi_rec,soi=as.matrix(soi),rec=as.matrix(rec))

tibble_soi_rec=tibble(df_soi_rec)
tsibble::duplicates(tibble_soi_rec, key = NULL, index=Fecha_soi_rec)   ##Mirar si hay registros duplicados

## # A tibble: 0 × 3
## # ℹ 3 variables: Fecha_soi_rec <date>, soi <dbl>, rec <dbl>

tibble_soi_rec%>%plot_acf_diagnostics(Fecha_soi_rec,soi,.ccf_vars = rec,.lags = 36)

En ambas gráficas podemos ver que se muestran periodicidades en las correlaciones que corresponden a valores separados por 12 unidades. También podemos ver que las observaciones con 12 meses o un año de diferencia están fuertemente correlacionadas positivamente, al igual que las observaciones en múltiplos como 24, 36, 48,\(\cdots\) Las observaciones separadas por seis meses están correlacionadas negativamente, lo que muestra que las excursiones positivas tienden a asociarse con las excursiones negativas a los seis meses ver Shumway2017 capítulo 1.

Índice AMI

###AMI Del los libros H. Kantz and T. Schreiber: Nonlinear Time series Analysis (Cambridge university press) H. Abarbanel: Analysis of observed chaotic data (Springer, 1996) y NONLINEAR TIME SERIES ANALYSIS HOLGER KANTZ AND THOMAS SCHREIBER. Cambrige University Press 2003.

Ahora utilizaremos los paquetes nonlinearTseries y tseriesChaos para computar el average mutual information(AMI) o La información mutua promedio (AMI, la cual mide cuánto nos dice una variable aleatoria sobre otra, el cual se define como:

\[I(X;Y)=\sum_{i}\sum_{j}p(x_i,y_j)\log_2(\frac{p(x_i,y_j)}{p(x_i)p(y_j)}).\] En el contexto del análisis de series de tiempo, AMI ayuda a cuantificar la cantidad de conocimiento obtenido sobre el valor de \(X_{t+d}\) al observar \(X_t\). Equivalentemente, el AMI es una medida de qué tanto el conocimiento de \(X\) reduce la incertidumbre acerca de \(Y\). Esto implica que \(I(X,Y)=0\) si y sólo si \(X\) y \(Y\) son variables aletorias independientes. I(X; Y ) describe la información que la medición \(X_t\) en el tiempo \(t\) aporta a la medición \(X_{t+d}\) en el tiempo \(t + d\). Si se elige d como el valor alrededor del primer mínimo del AMI, entonces \(Y{t}\) e \(y_{t+d}\) son parcialmente pero no totalmente independientes.

Vamos a simular una serie de la forma \[x_t=\frac{x_{t-1}}{x_{t-12}^2+1}+\epsilon_t\] y a trabajar con la serie de linces Candienses del paquete astsa lynx.

library(nonlinearTseries)

## 
## Attaching package: 'nonlinearTseries'

## The following object is masked from 'package:fabletools':
## 
##     estimate

## The following object is masked from 'package:grDevices':
## 
##     contourLines

library(tseriesChaos)
et=rnorm(1100,0,1)
xt=rep(0,1100)
for(t in 13:1100)
  {
  xt[t]=(xt[t-1]-1)/(xt[t-12]^2+1)+et[t]
  }
xtsimul=as.ts(xt[101:1100])
length(xtsimul)

## [1] 1000

plot(xtsimul)

acf(xtsimul)

par(mar = c(3,2,3,2))
astsa::lag1.plot(xtsimul, 12,corr = F) 

nonlinearTseries::mutualInformation(xtsimul,lag.max = 100,n.partitions = 50,units = "Bits",do.plot = TRUE) #c("Nats", "Bits", "Bans")

## $time.lag
##   [1]   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
##  [19]  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35
##  [37]  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53
##  [55]  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71
##  [73]  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89
##  [91]  90  91  92  93  94  95  96  97  98  99 100
## 
## $mutual.information
##   [1] 4.7211355 0.9086685 0.7868768 0.7146046 0.7234251 0.7393485 0.7257831
##   [8] 0.7204692 0.7530212 0.7388565 0.6987539 0.7426234 0.8363724 0.7918314
##  [15] 0.7495687 0.7384157 0.7227720 0.7306765 0.7185728 0.7490403 0.7888126
##  [22] 0.7335948 0.7676440 0.7241823 0.7728012 0.7690717 0.7278799 0.7436484
##  [29] 0.7654490 0.7561503 0.7324961 0.7820421 0.7413698 0.7871444 0.7796870
##  [36] 0.7356812 0.7411850 0.7501907 0.7512091 0.7691925 0.7466313 0.7425748
##  [43] 0.7152309 0.7567790 0.7658219 0.7516127 0.7371519 0.7288299 0.7622234
##  [50] 0.7844343 0.7584129 0.7388515 0.7564102 0.7690676 0.7334833 0.7605943
##  [57] 0.7808753 0.7568963 0.7571249 0.7135538 0.7925615 0.7584517 0.7768902
##  [64] 0.7298138 0.7462323 0.7553250 0.7728154 0.7850680 0.7464423 0.7543733
##  [71] 0.7465533 0.7813147 0.7620542 0.7528362 0.7663679 0.8111562 0.7237969
##  [78] 0.7386509 0.7523882 0.7398758 0.7844187 0.7326143 0.7380267 0.7455960
##  [85] 0.7691433 0.7309565 0.7479580 0.7314371 0.7583339 0.8090117 0.7704502
##  [92] 0.7519965 0.7847387 0.7292170 0.7291928 0.7169975 0.7426216 0.7759557
##  [99] 0.7898280 0.7688683 0.7700155
## 
## $units
## [1] "Bits"
## 
## $n.partitions
## [1] 50
## 
## attr(,"class")
## [1] "mutualInf"

pacf(xtsimul)

tseriesChaos::mutual(xtsimul, partitions = 50, lag.max = 100, plot=TRUE)

plot(astsa::Lynx)

acf(astsa::Lynx)

par(mar = c(3,2,3,2))
astsa::lag1.plot(astsa::Lynx, 12) 

tseriesChaos::mutual(astsa::Lynx, partitions = 50, lag.max = 10, plot=TRUE)

Detección de cíclos y estacionalidades

Podemos usar la función TSstudio::ts_heatmap para crear un mapa de calor. Este es un gráfico tridimensional, en donde en el eje y están los meses, y en el eje x están los años. Note en este caso que los meses de Diciembre,Enero, Febrero y Marzo presentan los valores mas oscuros, es decir los valores mas grandes a lo largo de los años, en contraste con los meses de Mayo a Septiembre que presentan colores mas claros. Este es un típico comportamiento de la presencia de un cíclo estacional en la serie. El flujo de color es horizontal.

TSstudio::USgas

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 2000 2510.5 2330.7 2050.6 1783.3 1632.9 1513.1 1525.6 1653.1 1475.0 1567.8
## 2001 2677.0 2309.5 2246.6 1807.2 1522.4 1444.4 1598.1 1669.2 1494.1 1649.1
## 2002 2487.6 2242.4 2258.4 1881.0 1611.5 1591.4 1748.4 1725.7 1542.2 1645.9
## 2003 2700.5 2500.3 2197.9 1743.5 1514.7 1368.4 1600.5 1651.6 1428.6 1553.2
## 2004 2675.8 2511.1 2100.9 1745.2 1573.0 1483.7 1584.9 1578.0 1482.2 1557.2
## 2005 2561.9 2243.0 2205.8 1724.9 1522.6 1534.1 1686.6 1695.1 1422.5 1428.2
## 2006 2165.3 2144.4 2126.4 1681.0 1526.3 1550.9 1758.7 1751.7 1462.1 1644.2
## 2007 2475.6 2567.0 2128.8 1810.1 1559.1 1555.2 1659.9 1896.1 1590.5 1627.8
## 2008 2734.0 2503.4 2278.2 1823.9 1576.4 1604.2 1708.6 1682.9 1460.9 1635.8
## 2009 2729.7 2332.5 2170.7 1741.3 1504.0 1527.8 1658.0 1736.5 1575.0 1666.5
## 2010 2809.8 2481.0 2142.9 1691.8 1617.3 1649.5 1825.8 1878.9 1637.5 1664.9
## 2011 2888.6 2452.4 2230.5 1825.0 1667.4 1657.3 1890.5 1891.8 1655.6 1744.5
## 2012 2756.2 2500.7 2127.8 1953.1 1873.8 1868.4 2069.8 2008.8 1807.2 1901.1
## 2013 2878.8 2567.2 2521.1 1967.5 1752.5 1742.9 1926.3 1927.4 1767.0 1866.8
## 2014 3204.1 2741.2 2557.9 1961.7 1810.2 1745.4 1881.0 1933.1 1809.3 1912.8
## 2015 3115.0 2925.2 2591.3 2007.9 1858.1 1899.9 2067.7 2052.7 1901.3 1987.3
## 2016 3091.7 2652.3 2356.3 2083.8 1965.8 2000.7 2186.6 2208.4 1947.8 1925.2
## 2017 2914.2 2340.6 2523.7 1932.5 1892.5 1910.9 2142.1 2094.3 1920.9 2032.0
## 2018 3335.0 2705.9 2792.6 2346.3 2050.9 2058.7 2344.6 2307.7 2151.5 2279.1
## 2019 3399.9 2999.2 2899.9 2201.1 2121.0 2115.2 2407.5 2437.2 2215.6 2472.3
##         Nov    Dec
## 2000 1908.5 2587.5
## 2001 1701.0 2120.2
## 2002 1913.6 2378.9
## 2003 1753.6 2263.7
## 2004 1782.8 2327.7
## 2005 1663.4 2326.4
## 2006 1765.4 2122.8
## 2007 1834.5 2399.2
## 2008 1868.9 2399.7
## 2009 1776.2 2491.9
## 2010 1973.3 2714.1
## 2011 2031.9 2541.9
## 2012 2167.8 2503.9
## 2013 2316.9 2920.8
## 2014 2357.5 2679.2
## 2015 2249.1 2588.2
## 2016 2159.4 2866.3
## 2017 2357.7 3084.5
## 2018 2709.9 2993.1
## 2019

 ts_plot(USgas,
           title = "US Monthly Natural Gas consumption",
           Ytitle = "Billion Cubic Feet",
           Xtitle = "Year",
           Xgrid = TRUE,
           Ygrid = TRUE)

TSstudio::ts_heatmap(USgas,title = "Mapa de Calor - Consumo de Gas Natural en EEUU")

Voy enfocarme en un periodo de la serie de desempleo.

 unemployment <- window(USUnRate, start = c(1990,1))
   ts_plot(unemployment,
           title = "US Monthly Unemployment Rate",
           Ytitle = "Unemployment Rate (%)",
           Xtitle = "Year",
           Xgrid = TRUE,
Ygrid = TRUE)

ts_info(unemployment)

##  The unemployment series is a ts object with 1 variable and 360 observations
##  Frequency: 12 
##  Start time: 1990 1 
##  End time: 2019 12

• El primer ciclo ocurrió entre 1990 y 2000, que estuvo cerca de un ciclo de 11 años.
• El segundo ciclo comenzó en 2000 y finalizó en 2007, que fue un ciclo de 7 años.
• Un tercer ciclo, que comenzó en 2007 y continúa a mayo de 2019 aún no se ha completado, lo que significa que se ha prolongado durante más de 12 años.

ts_heatmap(unemployment,title = "Mapa de Calor - Tasa de Desempleo EEUU Subserie")

#ts_heatmap(USUnRate,title = "Mapa de Calor - Tasa de Desempleo EEUU")

En este ejemplo, el flujo de color de la Tasa de Desempleo es vertical, lo que indica el estado del ciclo. En este caso, las franjas verticales más claras representan el final de un ciclo y el comienzo del siguiente. Asimismo, las franjas verticales más oscuras representan los picos del ciclo.

Vale la pena decir que es necesario eliminar la tendencia de la serie para pasar a detectar la estacionalidad.

Explorando mas herramientas para detección de Estacionalidad

Medidas descriptivas

USgas_df <- data.frame(year = floor(time(USgas)), month = cycle(USgas),USgas = as.numeric(USgas))

USgas_df$month <- factor(month.abb[USgas_df$month], levels = month.abb)

library(dplyr)
USgas_summary <- USgas_df %>%group_by(month) %>%summarise(mean= mean(USgas),sd = sd(USgas))
USgas_summary

## # A tibble: 12 × 3
##    month  mean    sd
##    <fct> <dbl> <dbl>
##  1 Jan   2806.  309.
##  2 Feb   2502.  223.
##  3 Mar   2325.  241.
##  4 Apr   1886.  175.
##  5 May   1708.  194.
##  6 Jun   1691.  216.
##  7 Jul   1864.  261.
##  8 Aug   1889.  237.
##  9 Sep   1687.  242.
## 10 Oct   1788.  260.
## 11 Nov   2015.  284.
## 12 Dec   2543.  278.

 library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

  plot_ly (data = USgas_summary, x = ~ month, y = ~ mean, type = "bar", name   = "Mean") %>%
   layout (title = "USgas - Monthly Average", yaxis =list(title = "Mean",   range = c(1500, 2700)))

  monthplot(USgas)

  ####También lo podemos hacer usando el paquete feasts y el objeto tsibble
  USgas_tsbl=as_tsibble(USgas)
  require(feasts)
  USgas_tsbl=as_tsibble(USgas)
  USgas_tsbl%>%gg_subseries(value)###Puede usar el argumento period=12 y da el mismo resultado, lo que significa es que se pueden agrupar las observaciones que están cada 12.

Note que basados en las estadísticas descriptivas podemos ver que las medias son distintas para algunos meses, incluso sin quitar la tendencia. A una misma conclusión llegamos basados en los gráficos. Estas son típicas características de que hay presente un ciclo estacional en la serie.

Explorando múltiples estacionalidades

Pueden haber múltiples estacionalidades(por lo general ocurren en series de alta frecuencia: diaria, cada hora, cada media y asi sucesivamente.) o una única estacionalidad, o inclusive cíclos que no se ven con facilidad.

Vamos a trabajar la base datos relacionada con el sistema nacional de transmisión de electricidad del Reino Unido. Más específicamente con la demanda de energía de forma horaria.

library(UKgrid)
require(TSstudio)
require(timetk)
require(feasts)
require(tsibble)
require(plotly)
UKgrid_xts <- extract_grid(type = "xts",
                              columns = "ND",
                              aggregate = "hourly",
                              na.rm = TRUE)
#extract_grid solo funciona para el conjunto de datos UKgrid
ts_plot(UKgrid_xts,
            title = "National Hourly Demand UK Grid",
            Ytitle = "Megawatts",
            Xtitle = "Year",
            Xgrid = TRUE,
            Ygrid = TRUE)

UKgrid_tstbl <- extract_grid(type = "tsibble",
                              columns = "ND",
                              aggregate = "hourly",
                              na.rm = TRUE)

UKgrid_tbl <-na.omit(as_tibble(UKgrid_tstbl))

Como indicamos anteriormente, la primera indicación de la posible existencia de múltiples patrones estacionales en la serie es una frecuencia alta, como por diaria, horaria y en minutos. En esos casos, hay más de una forma de establecer la frecuencia de la serie. Por ejemplo, si capturamos una serie de tiempo de frecuencia diaria, la frecuencia de la serie se puede configurar de la siguiente manera: * Diariamente (o 365), asumiendo que el ciclo más apropiado es un año completo. * Días de semana (o 7) siempre que la oscilación del día de la semana sea más dominante que la del ciclo de año completo.

Al utilizar estadísticas descriptivas con este tipo de series, tendrá sentido aplicar este método para cada frecuencia potencial de la serie (o al menos las principales) con el fin de examinar si existe una indicación del patrón estacional.

Por ejemplo, UKgrid es una serie de tiempo por horas, que la marca automáticamente como sospechosa de tener múltiples patrones estacionales. Potencialmente, como se mencionó anteriormente, la demanda horaria de electricidad podría tener tres patrones estacionales diferentes:

• Horaria: Este es probablemente el principal patrón estacional de la serie, ya que existe una relación directa entre la demanda de energía eléctrica y la hora del día (hay alta demanda durante el día y baja demanda durante la noche o al contrario).
• Día de la semana: La demanda de electricidad a lo largo del día se deriva, potencialmente, del día de la semana. Tendría sentido esperar un alto consumo durante los días laborables y una menor tasa de consumo durante el fin de semana.
• Mensual: como los patrones climáticos varían a lo largo del año, la cantidad de luz del día y otros factores estacionales podrían afectar la demanda de electricidad.

Usando el paquete lubridate crearemos las características.

library(xts)
UKgrid_df <- data.frame(time = zoo::index(UKgrid_xts), UKgrid=as.numeric(UKgrid_xts))
str(UKgrid_df)

## 'data.frame':    127296 obs. of  2 variables:
##  $ time  : POSIXct, format: "2005-04-01 00:00:00" "2005-04-01 01:00:00" ...
##  $ UKgrid: num  65080 68207 69172 66769 64660 ...

Ahora crearemos características estacionales basados en los periodos que deseamos explorar, por ejemplo hora del día,o día de la semana, o mes del año.

library(lubridate)
UKgrid_df$hour <- hour(UKgrid_df$time)
UKgrid_df$weekday <- wday(UKgrid_df$time, label = TRUE, abbr = TRUE)
UKgrid_df$month <- factor(month.abb[month(UKgrid_df$time)], levels =   month.abb)
head(UKgrid_df)

##                  time UKgrid hour weekday month
## 1 2005-04-01 00:00:00  65080    0     Fri   Apr
## 2 2005-04-01 01:00:00  68207    1     Fri   Apr
## 3 2005-04-01 02:00:00  69172    2     Fri   Apr
## 4 2005-04-01 03:00:00  66769    3     Fri   Apr
## 5 2005-04-01 04:00:00  64660    4     Fri   Apr
## 6 2005-04-01 05:00:00  65209    5     Fri   Apr

Vamos a empezar las exploraciones analizando el ciclo horario.

 UKgrid_hourly <- UKgrid_df %>%
    dplyr::group_by(hour) %>%
    dplyr::summarise(mean = mean(UKgrid, na.rm = TRUE), sd = sd(UKgrid, na.rm
= TRUE))
str(UKgrid_hourly)

## tibble [24 × 3] (S3: tbl_df/tbl/data.frame)
##  $ hour: int [1:24] 0 1 2 3 4 5 6 7 8 9 ...
##  $ mean: num [1:24] 55622 54343 52920 51546 50509 ...
##  $ sd  : num [1:24] 9259 9482 9485 9165 8557 ...

UKgrid_hourly

## # A tibble: 24 × 3
##     hour   mean     sd
##    <int>  <dbl>  <dbl>
##  1     0 55622.  9259.
##  2     1 54343.  9482.
##  3     2 52920.  9485.
##  4     3 51546.  9165.
##  5     4 50509.  8557.
##  6     5 51747.  8657.
##  7     6 59048. 10489.
##  8     7 68627. 13347.
##  9     8 73581. 13003.
## 10     9 76320. 12557.
## # ℹ 14 more rows

Vamos ahora a hacer la gráfica de la media y la desviación estándar con base en las horas. Note que esas gráficas en escalas diferentes, así que hay que usar un gráfico especial.

require(plotly)
 plot_ly(UKgrid_hourly) %>%
      add_lines(x = ~ hour, y = ~ mean, name = "Media") %>%
      add_lines(x = ~ hour, y = ~ sd, name = "Desviación Estándar", yaxis =
   "y2",
                line = list(color = "red", dash = "dash", width = 3)) %>%
      layout(
        title = "La demanda nacional de electricidad - Promedio horario vs. Desviación Estándar",
        yaxis = list(title = "Media"),
        yaxis2 = list(overlaying = "y",
                      side = "right",
                      title = "Desviación Estándar"
        ),
        xaxis = list(title="Hora del Día"),
        legend = list(x = 0.05, y = 0.9),
        margin = list(l = 50, r = 50)
)

Qué podemos destacar del gráfico y de las estadísticas descriptivas?

• Hay poca demanda durante la noche (entre la medianoche y las 6 a.m.) y una alta demanda entre las horas de la mañana y la tarde.

• Existe una fuerte correlación entre la demanda promedio y su desviación estándar.

• La relativamente baja desviación estándar de la demanda promedio durante la noche podría indicar que existe un fuerte efecto subestacional durante esas horas además de la estacionalidad horaria. Esto debería tener sentido, ya que son horas de sueño normales y, por lo tanto, en promedio, la demanda es razonablemente la misma durante los días de semana.

• Por otro lado, la alta desviación estándar a lo largo de las horas de alta demanda podría indicar que la demanda se distribuye de manera diferente en diferentes vistas de periodicidad (como día de la semana o mes del año).

Vamos a explorar el último punto, viendo la demanda a la madrugada(3 a.m) y empezando el día(9 a.m) con respecto al día de la semana.

   UKgrid_weekday <- UKgrid_df %>%
      dplyr::filter(hour == 3 | hour == 9) %>%
    dplyr::group_by(hour, weekday) %>%
    dplyr::summarise(mean = mean(UKgrid, na.rm = TRUE),
                     sd = sd(UKgrid, na.rm = TRUE))

## `summarise()` has grouped output by 'hour'. You can override using the
## `.groups` argument.

UKgrid_weekday$hour <- factor(UKgrid_weekday$hour)
    plot_ly(data = UKgrid_weekday, x = ~ weekday, y = ~ mean, type =
   "bar",color = ~ hour) %>%
      layout(title = "The Hourly Average Demand by Weekday",
             yaxis = list(title = "Mean", range = c(30000, 75000)),
             xaxis = list(title = "Weekday"))

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

En el gráfico de barras anterior podemos ver que la demanda de electricidad a las 3 a.m. es relativamente estable durante todos los días de la semana, con una ligera diferencia entre el promedio durante los días laborables y los días del fin de semana (alrededor de un 2% diferente). Por otro lado, existe una diferencia entre la demanda del día laborable y el fin de semana a las 9 a.m. (es decir, la demanda del lunes es en promedio un 28% superior a la del domingo). Como era de esperar, esos resultados se alinearon con nuestras expectativas anteriores.

Ahora podemos aprovechar esos conocimientos para examinar si existe un patrón estacional mensual en la serie. Ahora seleccionaremos las mismas horas (3 a.m. y 9 a.m.); sin embargo, esta vez agruparemos estos datos por mes (en lugar de días laborables):

 UKgrid_month <- UKgrid_df %>%
      dplyr::filter(hour == 3 | hour == 9) %>%
    dplyr::group_by(hour, month) %>%
    dplyr::summarise(mean = mean(UKgrid, na.rm = TRUE),
                     sd = sd(UKgrid, na.rm = TRUE))

## `summarise()` has grouped output by 'hour'. You can override using the
## `.groups` argument.

UKgrid_month$hour <- factor(UKgrid_month$hour)
    plot_ly(data = UKgrid_month, x = ~ month, y = ~ mean, type = "bar",color =
   ~ hour) %>%
      layout(title = "The Hourly Average Demand by Month",
             yaxis = list(title = "Mean", range = c(30000, 75000)),
             xaxis = list(title = "Month"))

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

Podemos ver en el gráfico de barras del resumen de agregación mensual que, en promedio, la demanda durante la noche (3 a.m.) y la mañana (9 a.m.) varía a lo largo de los meses del año. Además, hay un cambio significativo en la demanda durante la noche en comparación con la agregación entre semana. La variación de la serie de mes a mes indica la existencia de estacionalidad mensual en la serie.

Usando gráficos de densidades y Boxplot para explorar la estacionalidad.

Otro enfoque para analizar patrones estacionales en datos de series de tiempo es trazar la distribución de las unidades de frecuencia mediante el uso de histogramas o gráficos de densidad. Esto nos permitirá examinar si cada unidad de frecuencia tiene una distribución única que puede distinguirla del resto de unidades. Para esto usaremos el paquete ggplot2. Vamos empezar con la serie USgas.

require(timetk)
UKgrid_tbl

## # A tibble: 127,282 × 2
##    TIMESTAMP              ND
##    <dttm>              <int>
##  1 2005-04-01 00:00:00 65080
##  2 2005-04-01 01:00:00 68207
##  3 2005-04-01 02:00:00 69172
##  4 2005-04-01 03:00:00 66769
##  5 2005-04-01 04:00:00 64660
##  6 2005-04-01 05:00:00 65209
##  7 2005-04-01 06:00:00 72376
##  8 2005-04-01 07:00:00 82679
##  9 2005-04-01 08:00:00 88706
## 10 2005-04-01 09:00:00 91108
## # ℹ 127,272 more rows

#UKgrid_tbl%>%plot_seasonal_diagnostics(.date_var = TIMESTAMP,.value = ND,.feature_set = c("hour","week","wday.lbl","month.lbl"),.geom="boxplot") ##Chequear si se desea un diagrama de caja o violín .geom 

UKgrid_diff<-diff(UKgrid_tbl$ND)
which(is.na(UKgrid_diff))

## integer(0)

UKgrid_diff_sub=na.omit(UKgrid_diff)
str(UKgrid_diff_sub)

##  int [1:127281] 3127 965 -2403 -2109 549 7167 10303 6027 2402 209 ...

UKgrid_tbl%>%mutate(diff_ND=ND-dplyr::lag(ND))%>%plot_time_series(TIMESTAMP,diff_ND,.smooth=F)

UKgrid_tbl%>%mutate(diff_ND=ND-dplyr::lag(ND))%>%plot_seasonal_diagnostics(.date_var = TIMESTAMP,.value = diff_ND,.feature_set = c("hour"),.geom="boxplot")

## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).

UKgrid_tbl%>%mutate(diff_ND=ND-dplyr::lag(ND))%>%plot_seasonal_diagnostics(.date_var = TIMESTAMP,.value = diff_ND,.feature_set = c("wday.lbl"),.geom="boxplot")

## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).

UKgrid_tbl%>%mutate(diff_ND=ND-dplyr::lag(ND))%>%plot_seasonal_diagnostics(.date_var = TIMESTAMP,.value = diff_ND,.feature_set = c("month.lbl"),.geom="boxplot")

## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).

UKgrid_tstbl%>%gg_subseries(ND,period=24)

UKgrid_tstbl%>%gg_season(ND,period=24)

## Warning: Removed 14 rows containing missing values (`geom_line()`).

UKgrid_tstbl%>%gg_season(ND,period=24*7)

## Warning: Removed 14 rows containing missing values (`geom_line()`).

Usando Regresión para Descubrir un cíclo

Vamos a considerar la serie \[x_t =A\cos(2\pi\omega t+\varphi)+w_t,\] y una simulación de tamaño \(T=500\) para generar dos series, en donde las varianzas de los ruidos son diferentes en cada caso. Más específicamente \(\omega=\frac{1}{50}, A=2, \varphi=0.6\pi\) y \(\sigma_w=1,5.\)

cs = 2*cos(2*pi*1:500/50 + .6*pi);  w = rnorm(500,0,1)
par(mfrow=c(3,1), mar=c(3,2,2,1), cex.main=1.5)
plot.ts(cs, main=expression(2*cos(2*pi*t/50+.6*pi)))
plot.ts(cs+w, main=expression(2*cos(2*pi*t/50+.6*pi) + N(0,1)))
plot.ts(cs+5*w, main=expression(2*cos(2*pi*t/50+.6*pi) + N(0,25)))

Note que cuando la varianza del ruido es mas grande, es mas complicado ver el cíclo.

Vamos a re-escribir el proceso anterior usando la relación trigonométrica \(\cos(\alpha \pm \beta) = \cos(\alpha) \cos(\beta) \mp \sin(\alpha) \sin(\beta),\) con lo cual:

\[A \cos(2\pi\omega t + \varphi) = \beta_1 \cos(2\pi\omega t) + \beta_2 \sin(2\pi\omega t)\] donde \(\beta_1=A\cos(\varphi)\) y \(\beta_2=-A\sin(\varphi)\), así que el proceso queda escrito de la siguiente manera:

\[x_t = \beta_1 \cos(2\pi t/50) + \beta_2 \sin(2\pi t/50) + w_t.\] asimiendo que \(\omega=1/50\) es conocido. La transformación nos permitió convertir la regresión no lineal(en los parámetros) en una lineal. Así podemos estimar \(\beta_1\) y \(\beta_2\) usando una regresión lineal. Note también que si tomamos \(\beta_1^2+\beta_2^2\), entonces tenemos \(A^2=\beta_1^2+\beta_2^2\), así que, un estimador de la amplitud al cuadrado es \(\hat{A}^2=\hat{\beta}_1^2+\hat{\beta}_2^2.\) Note que los valores verdaderos de \(\beta_1=\) y \(\beta_2=\)

beta1teor=2*cos(.6*pi)
beta2teor=-2*sin(.6*pi)
set.seed(90210) # fijamos una semilla para que nos dé el mismo resultado
x = 2*cos(2*pi*1:500/50 + .6*pi) + rnorm(500,0,5)
z1 = cos(2*pi*1:500/50)
z2 = sin(2*pi*1:500/50)
summary(fit <- lm(x~0+z1+z2)) # sin intercepto.

## 
## Call:
## lm(formula = x ~ 0 + z1 + z2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8584  -3.8525  -0.3186   3.3487  15.5440 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## z1  -0.7442     0.3274  -2.273   0.0235 *  
## z2  -1.9949     0.3274  -6.093 2.23e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.177 on 498 degrees of freedom
## Multiple R-squared:  0.07827,    Adjusted R-squared:  0.07456 
## F-statistic: 21.14 on 2 and 498 DF,  p-value: 1.538e-09

beta1teor

## [1] -0.618034

beta2teor

## [1] -1.902113

par(mar = c(2,2,2,2))
par(mfrow=c(2,1))
plot.ts(x)
plot.ts(x, col=8, ylab=expression(hat(x)))
  lines(fitted(fit), col=2)

El periodograma

El periodograma es una herramienta que permite detectar el valor de la frecuencia \(\omega\) que en general es desconocida.

Exploración de Múltiples Cíclos y el Periodograma

Vamos a explorar una herremienta que sirve para describir funciones periódicas generales. Fourier mostró que toda función periódica puede representarse como una suma de funciones sinusoidales de distinta amplitud y frecuencia. La idea entonces es generalizar el análisis anterior para un ciclo a suma de funciones armónicas con distintas frecuencias.

Dada la serie de longitud \(T\), se denomina periodo básico o de Fourier, a las fracciones exactas del tamaño muestral, es decir \[s_j=\frac{T}{j},\ \ j=1,2,\cdots,T/2.\]

El valor máximo de periodo básico se obtiene para \(j=1\) y es \(T\), es decir, observamos la onda una sola vez. El mínimo se obtiene para \(j=T/2\) y es 2 ya que no se puede observar periodos que duren menos de 2 observaciones. Como antes, suele usarse las frecuencias en lugar de los periodos para el ajuste de los ciclos, así que las frecuencias básicas o de Fourier se define como \[f_j=\frac{j}{T}\ \ j=1,2,\cdots,T/2.\]

Así, \(1/T\leq f_j\leq 1/2\), y el valor máximo sería \(f_j=0.5\) y se conoce como la frecuencia Nyquist. Con esto, una serie de tiempo \(Z_t\) que es periódica se puede representar mediante

\[Z_t=\mu+\sum_{j=1}^{T/2}A_j\sin(\omega_j t)+ \sum_{j=1}^{T/2}B_j\cos(\omega_j t)+a_t.\]

Note que éste modelo tiene tantos parámetros como observaciones, con lo cual hay que buscar un procedimiento para seleccionar las frecuencias que deben ser incluidas para explicar la evolución de la serie. Para esto se utilizará la herramienta conocida como el periodograma. Se puede verificar que la contribución de una onda a la varianza es dada por \(\frac{\hat{R}^2}{2}\), así que ondas asociadas con amplitudes grandes serán importantes en la explicación de la variabilidad de la serie, mientras que ondas asociadas con amplitudes bajas serán poco importantes. De forma análoga a como se mostró en la sección anterior, los estimadores de los coeficientes \(A_j\) y \(B_j\) están dados por \[\hat{A}_j=\frac{2}{T}\sum_{t=1}^{T}Z_t\sin(\omega_j t)\]
\[\hat{B}_j=\frac{2}{T}\sum_{t=1}^{T}Z_t\cos(\omega_j t)\] y \[\hat{R}_j=\hat{A}_j^2+\hat{B}_j^2,\] donde \(\omega_j=2\pi f_j\) y \(f=\frac{j}{T}\), con lo cual \[TS_z^2=\frac{T}{2}\sum_{j=1}^{T/2}\hat{R}_j^2.\]

Se le da el nombre de periodograma a la representación de cada \(\frac{T\hat{R}_j^2}{2}\) en función de la frecuencia \(\omega_j\) o \(f_j\). Así \[\begin{equation} \label{periodograma} I(f_j)=\frac{T\hat{R}_j^2}{2},\ \ \text{con}\ 1/T\leq f_j\leq 1/2 \end{equation}\] y \[\bar{I}=\sum_{j=1}^{T/2}\hat{R}_j^2=S_z^2.\]

Nos enfocaremos en las frecuencias básicas únicamente, lo cual no es restrictivo si \(T\) es muy grande, puesto que el número de frecuencias básicas será muy grande y siempre existirá alguna frecuencia básica muy próxima a la que puede interesarnos. El periodograma puede verse como una herramienta para la detección de posibles ciclos(ocultos) deterministas en una serie temporal. Por ejemplo, en una serie mensual estacional de periodo \(s=12\), esperamos encontrar un valor alto de periodograma para \(f=1/12\), pero también para \(f=j/12\), es decir, \(1/6,1/4,1/3\) que son armónicos del periodo estacional. Por otro lado, la serie puede tener otros ciclos no necesariamente ligado al periodo estacional, siendo el periodograma una buena herramienta para detectar estos posibles componentes. \

Podemos considerar la amplitud calculada para la frecuencia \(f_j\) como un promedio de las amplitudes existentes en las frecuencias situadas en el intervalo \(f_j\pm \frac{1}{2}T\). Con esto, se puede obtener el periodograma suavizado construyendo rectángulos con centro \(f_j\), base igual a \(1/T\) y alturas \(I(f_j)=T\hat{R}_j/2\). Otra forma de suavizar el periodograma es \[I(f)=\sum_{f_i-q}^{fi+q}p_iI(f_i)\ \ 0\leq f \leq 0.5,\]

donde \(q\) representa la ventana usada y los \(p_i\) son los pesos simétricos y no modifica el valor de la varianza \(S^2_z\) que es área bajo la curva.

Vamos a retomar el ejemplo anterior donde simulamos una serie de tamaño \(T=500\) con frecuencia \(\omega=1/50\).

plot.ts(x)

spectrum(x,log='no')
abline(v=1/50, lty=2,col="red")

###Otra forma de computarlo es vía la transformada de Fourier

Px = Mod(fft(x))^2
Freq=0:499/500
plot(Freq, Px, type='l')

u = which.max(Px[1:250])
sprintf("El valor de la frecuencia donde se máximiza el periodograma es %s",Freq[u])

## [1] "El valor de la frecuencia donde se máximiza el periodograma es 0.02"

Note que en cualquiera de las dos formas de obtener el periodograma, hemos descubierto la frecuencia \(\omega=1/50\) como frecuencia príncipal.

En ocasiones es necesario suavizar el periodograma porque encontramos varios picos que en verdad no son valores grandes.

spectrum(x,log='no',span=5)

#span vector of odd integers giving the widths of modified Daniell smoothers to be used to smooth the periodogram.

spectrum(x,log='no',span=c(5,5))

spectrum(x,log='no',span=c(2,2))

Qué sucede si no consideramos un ciclo determinístico a través de funciones sinusoidales sino que lo consideramos a través de variables Dummy? Es decir por ejemplo una componente estacional de periodo 4, modelado a través de variables dummy. \[Y_t=\delta_1\gamma_{1,t}+\delta_2\gamma_{2,t}+\delta_3\gamma_{3,t}+\delta_4\gamma_{4,t}+\epsilon_{t}\] donde las variables dummy \(\gamma_{i,t}, i=1,2,3,4\) indican si se está en el trimestre \(i\). Puede en periodograma detactar este tipo de estacionalidades obtenidas a través de dummy estacionales?

library(dplyr)
library(rstanarm)

## Loading required package: Rcpp

## This is rstanarm version 2.21.4

## - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!

## - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.

## - For execution on a local, multicore CPU with excess RAM we recommend calling

##   options(mc.cores = parallel::detectCores())

## 
## Attaching package: 'rstanarm'

## The following object is masked from 'package:car':
## 
##     logit

## The following objects are masked from 'package:VGAM':
## 
##     cauchy, dirichlet, exponential, laplace, logit

library(rstan)

## Loading required package: StanHeaders

## rstan (Version 2.21.8, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)

## 
## Attaching package: 'rstan'

## The following object is masked from 'package:tidyr':
## 
##     extract

library(ggplot2)
library(bayesplot)

## This is bayesplot version 1.10.0

## - Online documentation and vignettes at mc-stan.org/bayesplot

## - bayesplot theme set to bayesplot::theme_default()

##    * Does _not_ affect other ggplot2 plots

##    * See ?bayesplot_theme_set for details on theme setting

s=4
deltas=seq(1:s)

####Ejemplo Fácil Dummy+ruido

T=200*s
epsilon=rnorm(T,0,1)
l=diag(s)
X=matrix(rep(t(l),T/s),ncol=ncol(l),byrow=TRUE)
Y=X%*%deltas+epsilon
plot(as.ts(Y))

acf(as.ts(Y))

#simul=data.frame(resp=Y,diseno=X)
#salida_dummy_ruido=stan_glm(resp~. -1,data=simul)
PeriodgramaY=spectrum(Y,log='no')

ubicacionY=which.max(PeriodgramaY$spec)
sprintf("El valor de la frecuencia donde se máximiza el periodograma para la serie es: %s",PeriodgramaY$freq[ubicacionY])

## [1] "El valor de la frecuencia donde se máximiza el periodograma para la serie es: 0.25"

sprintf("El periodo correspondiente es aporximadamente: %s",1/PeriodgramaY$freq[ubicacionY])

## [1] "El periodo correspondiente es aporximadamente: 4"

Veamos ahora con una serie de tiempo real. Usaremos la serie SOi y Recruitment para ilustrar su uso. Recordemos que ambas series son mensuales, así que las frecuencias que serán dibujadas en el eje x vendrán en múltiplos de \(\frac{1}{12}\). Usaremos ahora la función mvspec del paquete astsa.

library(astsa)
data(soi)
plot(soi)

soi.per = astsa::mvspec(soi, log="no")

aux.per=spectrum(soi[1:453])
ubicacionsoi.aux=which.max(aux.per$spec)
ubicacionsoi=which.max(soi.per$spec)
sprintf("El valor de la frecuencia donde se máximiza el periodograma para SOI es: %s",aux.per$freq[ubicacionsoi.aux])

## [1] "El valor de la frecuencia donde se máximiza el periodograma para SOI es: 0.0833333333333333"

sprintf("El valor de la frecuencia donde se máximiza el periodograma para SOI es: %s",soi.per$freq[ubicacionsoi])

## [1] "El valor de la frecuencia donde se máximiza el periodograma para SOI es: 1"

n_soi <- length(soi.per$spec)
valor_seg_per=sort(soi.per$spec,partial=n_soi-1)[n_soi-1]
ubica_segundo_soi=which(soi.per$spec==valor_seg_per)

sprintf("El valor de la frecuencia donde se alcanza el segundo máximo para el periodograma para SOI es: %s",soi.per$freq[ubica_segundo_soi])

## [1] "El valor de la frecuencia donde se alcanza el segundo máximo para el periodograma para SOI es: 0.2"

abline(v=soi.per$freq[ubicacionsoi], lty=2,col="blue")
abline(v=soi.per$freq[ubica_segundo_soi], lty=2,col="blue")

Note que la primera frecuencia mas alta se encuentra en \(\omega=1\cdot \frac{1}{12}=\frac{1}{12}\), el cual corresponde al obvio ciclo anual (periodo de 12 meses). La segunda corresponde a \(\omega=\frac{1}{5} \cdot \frac{1}{12}=\frac{1}{60}\), el cual corresponde a un posible ciclo de periodo \(1/(1/60)= 60\) meses, es decir 5 años, lo cual puede ser debido al fenómeno de El Niño.

data("rec")
plot(rec)

rec.per = astsa::mvspec(rec, log="no")
ubicacionrec=which.max(rec.per$spec)
sprintf("El valor de la frecuencia donde se máximiza el periodograma para REC es: %s",rec.per$freq[ubicacionrec])

## [1] "El valor de la frecuencia donde se máximiza el periodograma para REC es: 1"

n_rec <- length(rec.per$spec)
valor_seg_per_rec=sort(rec.per$spec,partial=n_rec-1)[n_rec-1]
ubica_segundo_rec=which(rec.per$spec==valor_seg_per_rec)

sprintf("El valor de la frecuencia donde se alcanza el segundo máximo para el periodograma para REC es: %s",rec.per$freq[ubica_segundo_rec])

## [1] "El valor de la frecuencia donde se alcanza el segundo máximo para el periodograma para REC es: 0.25"

abline(v=rec.per$freq[ubicacionrec], lty=2,col="blue")
abline(v=rec.per$freq[ubica_segundo_rec], lty=2,col="blue")

Note que la primera frecuencia mas alta se encuentra en \(\omega=1\cdot \frac{1}{12}=\frac{1}{12}\), el cual corresponde al obvio ciclo anual (periodo de 12 meses), igual que para el casi del indice soi. La segunda corresponde a \(\omega=\frac{1}{4} \cdot \frac{1}{12}=\frac{1}{48}\), el cual corresponde a un posible ciclo de periodo \(1/(1/48)= 48\) meses, es decir 4 años. lo cual puede ser debido al fenómeno de El Niño. Esta actividad de banda ancha sugiere que el posible ciclo de El Niño es irregular, pero tiende a rondar los cuatro años en promedio.

Tarea : Hacer lo mismo pero ahora usar la serie de pasajeros diferenciada, al igual que la serie UKgrid, mas específicamente la serie UKgrid_xts y la serie taylor del paquete forecast.

Componente Estacional Estocástica

La componente estacional \(S_t^{(s)}\) puede ser estocástica y no determinística. Consideremos el primer caso, donde la evolución es estacionaria, es decir los estacionales no son constante, pero siguen un proceso estacionario, oscilando alrededor de un valor medio de acuerdo con la representación propuesta en [@Pena2010]: \[S_t^{(s)}=\mu^{(s)}+\nu_t\] donde \(\nu_t\) es un proceso estacionario de media cero y \(\mu^{(s)}\) es una constante que depende del mes o de las componentes de los \(s\) periodos.

Otra forma de introducir la estacionalidad de forma estocástica es permitir que sea cambiante en el tiempo sin ningún valor fijo medio. En este caso la estacionalidad puede proponerse que siga un proceso no estacionario por ejemplo que siga una caminata aleatoria: \[S_t^{(s)}=S_{t-s}^{(s)}+\nu_t\] siendo \(\nu_t\) es un proceso estacionario de media cero.

library(sarima) 

## 
## Attaching package: 'sarima'

## The following object is masked from 'package:rstanarm':
## 
##     se

## The following object is masked from 'package:astsa':
## 
##     sarima

## The following object is masked from 'package:stats':
## 
##     spectrum

x <- ts(sarima::sim_sarima(n=144, model = list(siorder=1,
                                      nseasons=12, sigma2 = 1)),frequency = 12 )#(1-B)^{12} X_t = e_t
plot(x)

acf(x,lag.max = 48)

monthplot(x)

ts_seasonal(x, type = "all")

spectrum(x,log='no',span=5)

## Estimated spectral density
##     series: x
##     method: Smoothed Periodogram
##     nseasons: 12
##     frequency range: (-6,6]
## 
##     sort method for the table: decreasing magnitudes
## 
##               freq      spec    % Total     Cum. %
##     [1,] 1.0000000 0.5488446 0.09121551 0.09121551
##     [2,] 0.9166667 0.5259173 0.08740510 0.17862061
##     [3,] 1.0833333 0.5008903 0.08324572 0.26186633
##     [4,] 0.8333333 0.3152458 0.05239244 0.31425877
##     [5,] 2.0833333 0.3137395 0.05214210 0.36640087
##     [6,] 2.0000000 0.3113759 0.05174928 0.41815014

astsa::mvspec(x, log="no")

library(sarima) 
x <- ts(sim_sarima(n=288,model=list(sar=0.8, nseasons=12, sigma2 = 1)), frequency = 12 )#SAR(1)
plot(x)

acf(x)

monthplot(x)

ts_seasonal(x, type = "all")

spectrum(x,log='no',span=5)

## Estimated spectral density
##     series: x
##     method: Smoothed Periodogram
##     nseasons: 12
##     frequency range: (-6,6]
## 
##     sort method for the table: decreasing magnitudes
## 
##              freq      spec    % Total     Cum. %
##     [1,] 2.000000 0.7349725 0.06086384 0.06086384
##     [2,] 2.041667 0.6623596 0.05485069 0.11571453
##     [3,] 1.958333 0.6055224 0.05014394 0.16585846
##     [4,] 2.083333 0.4523107 0.03745632 0.20331478
##     [5,] 5.000000 0.3843284 0.03182663 0.23514141
##     [6,] 1.916667 0.3716113 0.03077352 0.26591493

astsa::mvspec(x, log="no")

raices=polyroot(c(1,rep(0,11),-0.3))
raices

##  [1]  0.5527684+0.9574230i -0.9574230+0.5527684i -0.5527684-0.9574230i
##  [4]  0.9574230-0.5527684i  0.0000000+1.1055369i -1.1055369-0.0000000i
##  [7]  0.0000000-1.1055369i  1.1055369+0.0000000i -0.5527684+0.9574230i
## [10] -0.9574230-0.5527684i  0.5527684-0.9574230i  0.9574230+0.5527684i

conj_raices=Conj(raices)
raices[12]*conj_raices[12]

## [1] 1.222212+0i

Desestacionalizar o eliminación de la componente estacional.

Al igual como se hizo con la tendencia, hay al menos dos forma de eliminar la componente estacional: * Estimando la componente estacional y luego proceder quitarla. * Usando diferencias estacionales \((1-B^s)^D\), donde \(s\) es el periodo del ciclo estacional.

De las series con estacionalidad, usar el operador diferencia estacional y proceder a eliminar la componente estacional.De las series con estacionalidad, usar el operador diferencia estacional y proceder a eliminar la componente estacional usando la función diff(series,lag=s,differences=D).

spectrum(UKgrid_tbl$ND,log='no')

## Estimated spectral density
##     series: UKgrid_tbl$ND
##     method: Raw Periodogram
##     nseasons: 1
##     frequency range: (0,0]
## 
##     sort method for the table: decreasing magnitudes
## 
##                  freq      spec    % Total    Cum. %
##     [1,] 0.0416718750 23877.989 0.37309651 0.3730965
##     [2,] 0.0001171875  7230.329 0.11297478 0.4860713
##     [3,] 0.0000078125  6033.596 0.09427569 0.5803470
##     [4,] 0.0833437500  4327.080 0.06761115 0.6479581
##     [5,] 0.0059531250  3483.174 0.05442502 0.7023831
##     [6,] 0.0001093750  3032.342 0.04738072 0.7497639

spectrum(UKgrid_diff_sub,log='no')

## Estimated spectral density
##     series: UKgrid_diff_sub
##     method: Raw Periodogram
##     nseasons: 1
##     frequency range: (0,0]
## 
##     sort method for the table: decreasing magnitudes
## 
##                freq      spec    % Total    Cum. %
##     [1,] 0.04167188 22234.738 0.34742049 0.3474205
##     [2,] 0.08334375 15843.855 0.24756217 0.5949827
##     [3,] 0.12501563  3399.641 0.05311980 0.6481025
##     [4,] 0.20835938  2566.200 0.04009719 0.6881997
##     [5,] 0.16668750  2470.996 0.03860961 0.7268093
##     [6,] 0.12489844   672.260 0.01050415 0.7373134

PeriodgramaUKgrid_diff=spectrum(UKgrid_diff_sub,log='no')
ubicacionUKgrid_diff=which.max(PeriodgramaUKgrid_diff$spec)
sprintf("El valor de la frecuencia donde se máximiza el periodograma para la serie es: %s",PeriodgramaUKgrid_diff$freq[ubicacionUKgrid_diff])

## [1] "El valor de la frecuencia donde se máximiza el periodograma para la serie es: 0.041671875"

sprintf("El periodo correspondiente es aporximadamente: %s",1/PeriodgramaUKgrid_diff$freq[ubicacionUKgrid_diff])

## [1] "El periodo correspondiente es aporximadamente: 23.9970003749531"

###Ubicando el segundo valor mas grande
n.spec=length(PeriodgramaUKgrid_diff$spec)
valor_seg_per=sort(PeriodgramaUKgrid_diff$spec,partial=n.spec-1)[n.spec-1]
ubica_segundo_frec=which(PeriodgramaUKgrid_diff$spec==valor_seg_per)
1/PeriodgramaUKgrid_diff$freq[ubica_segundo_frec]

## [1] 11.9985

spectrum(UKgrid_tbl$ND,log='no')

## Estimated spectral density
##     series: UKgrid_tbl$ND
##     method: Raw Periodogram
##     nseasons: 1
##     frequency range: (0,0]
## 
##     sort method for the table: decreasing magnitudes
## 
##                  freq      spec    % Total    Cum. %
##     [1,] 0.0416718750 23877.989 0.37309651 0.3730965
##     [2,] 0.0001171875  7230.329 0.11297478 0.4860713
##     [3,] 0.0000078125  6033.596 0.09427569 0.5803470
##     [4,] 0.0833437500  4327.080 0.06761115 0.6479581
##     [5,] 0.0059531250  3483.174 0.05442502 0.7023831
##     [6,] 0.0001093750  3032.342 0.04738072 0.7497639

PeriodgramaUKgrid=spectrum(UKgrid_tbl$ND,log='no')
ubicacionUKgrid=which.max(PeriodgramaUKgrid$spec)
sprintf("El valor de la frecuencia donde se máximiza el periodograma para la serie es: %s",PeriodgramaUKgrid_diff$freq[ubicacionUKgrid_diff])

## [1] "El valor de la frecuencia donde se máximiza el periodograma para la serie es: 0.041671875"

sprintf("El periodo correspondiente es aporximadamente: %s",1/PeriodgramaUKgrid_diff$freq[ubicacionUKgrid])

## [1] "El periodo correspondiente es aporximadamente: 23.9970003749531"

###Ubicando el segundo valor mas grande
n.spec=length(PeriodgramaUKgrid$spec)
valor_seg_per=sort(PeriodgramaUKgrid$spec,partial=n.spec-1)[n.spec-1]
ubica_segundo_frec=which(PeriodgramaUKgrid$spec==valor_seg_per)

sprintf("El periodo correspondiente es aporximadamente: %s",1/PeriodgramaUKgrid$freq[ubica_segundo_frec])

## [1] "El periodo correspondiente es aporximadamente: 8533.33333333333"

#Está bastante cerca de 24*365=8760 el cual es un periodo anual

Tarea: Por favor hacer los análisis descritos anteriormente,es decir los gráficos de barras de las medias por horas, días de la semana, y meses para detectar la estacionalidad tanto para serie diferenciada, como para la serie cuando se estima la tendencia y luego se elimina(Use fitro de promedios móviles o STL).

 library(ggplot2)
    ggplot(USgas_df, aes(x = USgas)) +
      geom_density(aes(fill = month)) +
      ggtitle("USgas - Kernel Density Estimates by Month") +
      facet_grid(rows = vars(as.factor(month)))

Podemos ver algunos indicios de un patrón estacional en la serie, ya que las gráficas de densidad no se superponen entre sí (con la excepción de algunos meses consecutivos, como mayo y junio). Además, podemos ver que, durante algunos meses, la forma de las distribuciones es más plana con colas largas (principalmente durante los meses de invierno, noviembre, diciembre y enero). Sin embargo, no olvidemos el efecto de la tendencia o el crecimiento de un año a otro (como sabemos del capítulo anterior, la serie de gas de Estados Unidos tuvo una tendencia lineal desde el año 2010) ya que no la eliminamos de la serie. Repitamos este proceso; esta vez quitando la tendencia de la serie USgas antes de graficarla. Vamos a eliminarle la tendencia, posteriormente los explicaremos.

USgas_df$USgas_detrend <- USgas_df$USgas - decompose(USgas)$trend
    ggplot(USgas_df, aes(x = USgas_detrend)) +
      geom_density(aes(fill = month)) +
      ggtitle("USgas - Estimación de la densidad vía Kernel por mes") +
      facet_grid(rows = vars(as.factor(month)))

## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

## Warning: Removed 12 rows containing non-finite values (`stat_density()`).

Podemos ver que hay un comportamiento similar pero ahora las colas de las densidades estimadas son mas cortas.Pueden hacer lo mismo pero en vez de quitar la tendencia, hacemos una diferenciación.

En el caso de que la distribución de la mayoría de las unidades de frecuencia sea plana con una cola larga, podría ser una indicación de múltiples patrones estacionales en la serie. Regresemos a la serie UKgrid y tracemos las gráficas de densidad de 24 horas:

UKgrid_df$hour <- as.factor(UKgrid_df$hour)
    ggplot(UKgrid_df, aes(x = UKgrid)) +
      geom_density(aes(fill = hour)) +
      ggtitle("UKgrid - Kernel Density Estimates by Hour of the day") +
      facet_grid(rows = vars(as.factor(hour)))

## Warning: Removed 14 rows containing non-finite values (`stat_density()`).

Como observamos anteriormente con las tablas de resúmenes estadísticos, la distribución de la demanda neta de electricidad durante la noche es relativamente estable (de ahí la distribución no plana con colas cortas en contraposición a la distribución plana con cola larga durante el día). Si ahora hacemos un subconjunto con una de las horas durante el día y trazamos su distribución por el día de la semana, deberíamos esperar una superposición durante la noche y poder distinguir entre la distribución durante los días de la semana y el fin de semana, en contraposición a solo el día de la semana.

Por ejemplo, la siguiente gráfica representa la distribución de la demanda a las 9 a.m. a lo largo de los días de la semana. Puede ver que la distribución durante los días de la semana se distingue de la del fin de semana:

UKgrid_df$weekday <- as.factor(UKgrid_df$weekday)
    UKgrid_df %>% dplyr::filter(hour == 0) %>%
    ggplot(aes(x = UKgrid)) +
      geom_density(aes(fill = as.factor(weekday))) +
      ggtitle("UKgrid - Kernel Density Estimates by Hour of the day") +
      facet_grid(rows = vars(as.factor(weekday)))

Más herramientas del análisis de la estacionalidad

Vamos a usar el paquete forecast. Note que extraemos las subseries de los años. Lo cual demuestra un fuerte patrón estacional mensual.

library(forecast)
   ggseasonplot(USgas,year.labels=TRUE,continuous=TRUE)

ggseasonplot(USgas,  polar = TRUE)

Note que con el paquete TSstudio se puede hacer algo análogo.

ts_seasonal(USgas,type ="normal")

ts_seasonal(USgas, type = "cycle")

ts_seasonal(USgas, type = "box")

ts_seasonal(USgas, type = "all")

ts_seasonal con type=“cycle” añade un orden cronológico, en este caso por mes. Esto puede permitir la identificación de un patrón estacional sin tener que quitar la tendencia.

Cuando en el argumento type=“box”, elabora un gráfico de caja por unidad de frecuencia.

Finalmente cuando en el argumento type=“all”, elabora todos los gráficos anteriores.

El mapa de calor ya lo vimos anteriormente.

ts_heatmap(USgas, color = "Reds")

Gráficos basados en cuantiles también son de utilidad.De forma predeterminada, la función devuelve un gráfico de cuantiles de las unidades de frecuencia de la serie, donde la línea media representa la mediana y las líneas inferior y superior representan los percentiles 25 y 75.

ts_quantile(UKgrid)

## Warning in ts_quantile(UKgrid): The input object is a multiple time series
## object, by defualt will use only the first series as an input

El argumento period permite examinar si los patrones estacionales de la serie están cambiando cuando se usa un subconjunto de tiempo diferente. Esto le permite examinar si la serie tiene patrones estacionales adicionales. Por ejemplo, podemos trazar el ciclo de 24 horas de la serie UKgrid por día de la semana estableciendo el argumento de período en días de la semana:

ts_quantile(UKgrid, period = "weekdays", n = 2)

## Warning in ts_quantile(UKgrid, period = "weekdays", n = 2): The input object is
## a multiple time series object, by defualt will use only the first series as an
## input

Como vimos anteriormente con las gráficas de densidad, la demanda de electricidad durante el día es relativamente mayor durante los días de semana en comparación con los fines de semana. De la misma manera, puede trazar el ciclo de 24 horas por mes:

 ts_quantile(UKgrid, period = "monthly", n = 2)

## Warning in ts_quantile(UKgrid, period = "monthly", n = 2): The input object is
## a multiple time series object, by defualt will use only the first series as an
## input

La principal ventaja de las gráficas de cuantiles de múltiples períodos (por ejemplo, días de la semana, meses, etc.) sobre la gráfica de densidad que usamos anteriormente es que la primera representa todas las unidades de frecuencia (y, en el caso de UKgrid serie, con una frecuencia de 24 horas), mientras que el segundo representa una sola unidad de frecuencia (por ejemplo, la densidad de las 9 am durante los días de semana)

library(timetk)
library(broom)
UKgrid_tbl_diff<-UKgrid_tbl|>mutate(diff_ND=ND-dplyr::lag(ND))
num_harmonics=2
frequency=1/24



UKgrid_tbl_diff_Fourier<-UKgrid_tbl_diff|>mutate(
  C1_24=fourier_vec(TIMESTAMP,period=1/frequency,K=1,type="cos"),
  S1_24=fourier_vec(TIMESTAMP,period=1/frequency,K=1,type="sin")
  )

UKgrid_tbl_diff_Fourier2<-UKgrid_tbl_diff|>mutate(
  C1_24=fourier_vec(TIMESTAMP,period=1/frequency,K=1,type="cos"),
  S1_24=fourier_vec(TIMESTAMP,period=1/frequency,K=1,type="sin"),
  C2_24=fourier_vec(TIMESTAMP,period=1/frequency,K=2,type="cos"),
  S2_24=fourier_vec(TIMESTAMP,period=1/frequency,K=2,type="sin")
  )

####Un Armónico
Fourier1<-lm(diff_ND~C1_24+S1_24,data = UKgrid_tbl_diff_Fourier)

summary(Fourier1)

## 
## Call:
## lm(formula = diff_ND ~ C1_24 + S1_24, data = UKgrid_tbl_diff_Fourier)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7913.8 -2541.1  -644.6  1718.7 14207.4 
## 
## Coefficients:
##               Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)    -0.5401     9.6724   -0.056    0.955    
## C1_24       -2789.1503    13.6795 -203.892   <2e-16 ***
## S1_24        2377.7692    13.6781  173.838   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3451 on 127278 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.3606, Adjusted R-squared:  0.3606 
## F-statistic: 3.589e+04 on 2 and 127278 DF,  p-value: < 2.2e-16

plot(UKgrid_tbl_diff_Fourier$TIMESTAMP[2:49],UKgrid_tbl_diff_Fourier$diff_ND[2:49],type="l") 

lines(UKgrid_tbl_diff_Fourier$TIMESTAMP[2:49],Fourier1$fitted.values[1:48],type="l",col = 3) 

####Dos Armónicos Fourier
Fourier2<-lm(diff_ND~C1_24+S1_24+C2_24 +S2_24,data = UKgrid_tbl_diff_Fourier2)

summary(Fourier2)

## 
## Call:
## lm(formula = diff_ND ~ C1_24 + S1_24 + C2_24 + S2_24, data = UKgrid_tbl_diff_Fourier2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7917.1 -1802.0  -422.5  1429.7 11476.1 
## 
## Coefficients:
##               Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)    -0.8639     7.6550   -0.113     0.91    
## C1_24       -2789.7775    10.8264 -257.683   <2e-16 ***
## S1_24        2377.9248    10.8253  219.664   <2e-16 ***
## C2_24       -2965.0975    10.8262 -273.882   <2e-16 ***
## S2_24         327.6164    10.8255   30.263   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2731 on 127276 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5995, Adjusted R-squared:  0.5995 
## F-statistic: 4.763e+04 on 4 and 127276 DF,  p-value: < 2.2e-16

plot(UKgrid_tbl_diff_Fourier2$TIMESTAMP[2:49],UKgrid_tbl_diff_Fourier2$diff_ND[2:49],type="l") 

lines(UKgrid_tbl_diff_Fourier2$TIMESTAMP[2:49],Fourier2$fitted.values[1:48],type="l",col = 3) 

Suavizamiento en el contexto de Series de Tiempo

Promedio Móvil

El promedio móvil es un método es útil para descubrir ciertos rasgos en una serie de tiempo, como tendencias a largo plazo y componentes estacionales. En particular, si \(x_t\) representa las observaciones, entonces una forma de predecir predecir o estimar la tendencia de la serie es:

\[m_t=\sum_{j=-k}^{k}a_jx_{t-j},\] donde si\(a_j=a_{-j}\ge0\) y \(\sum_{j=-k}^{k}a_j=1\) se conoce como el promedio móvil simétrico de los datos.

# Filtro Boxcar
wgts = c(.5, rep(1,11), .5)/12   ###Los pesos de los filtros
soif = stats::filter(soi, sides=2, filter=wgts)
plot(soi)
lines(soif, lwd=2, col=4)

###Forma del Filtro
nwgts = c(rep(0,20), wgts, rep(0,20))
plot(nwgts, type="l", ylim = c(-.02,.1), xaxt='n', yaxt='n', ann=FALSE)

Suavizamiento Kernel

El suavizamiento kernel es un suavizador de promedio móvil que utiliza una función de ponderación, o kernel, para promediar las observaciones.Veamos ahora como queda el promedio móvil:

\[m_t=\sum_{i=1}^{n}w_i(t)x_i\] donde \[w_i(t)=K(\frac{t-i}{b})/\sum_{j=1}^{n}K(\frac{t-j}{b})\] son los pesos y \(K(\cdot)\) es una función kernel. Este estimador es llamado el Estimador de Nadaraya-Watson. Usaremos la función ksmooth

plot(soi)
lines(ksmooth(time(soi), soi, "normal", bandwidth=1), lwd=2, col=4)

#par(fig = c(.65, 1, .65, 1), new = TRUE) # the insert
gauss = function(x) { 1/sqrt(2*pi) * exp(-(x^2)/2) }
x = seq(from = -3, to = 3, by = 0.001)
plot(x, gauss(x), type ="l", ylim=c(-.02,.45), xaxt='n', yaxt='n', ann=FALSE)

Lowess

Otro enfoque para suavizar un gráfico de tiempo es la regresión del vecino más cercano. La técnica se basa en la regresión de k vecinos más cercanos, en la que uno usa solo los datos \(\{x_{t−k/2}, ..., x_t, ..., x_{t+k/2}\}\) para predecir \(x_t\) mediante regresión, y luego establece \(m_t = \hat{x}_t\). Primero, una cierta proporción de vecinos más cercanos a \(x_t\) se incluyen en un esquema de ponderación; los valores más cercanos a \(x_t\) en el tiempo obtienen más peso. Luego, se utiliza una regresión ponderada robusta para predecir \(x_t\) y obtener los valores suavizados \(mt\). Cuanto mayor sea la fracción de vecinos más cercanos incluidos, más suave será el ajuste. El R se usa la función lowess.

plot(soi)
lines(lowess(soi, f=.05), lwd=2, col=4)  # Ciclo El Niño 
lines(lowess(soi), lty=2, lwd=2, col=2)  # tendencia (con span por defecto)

Suavizamiento Splines

Una forma obvia de suavizar los datos sería ajustar una regresión polinomial en términos del tiempo. Por ejemplo, un polinomio cúbico tendría \(x_t = m_t + w_t\) donde \(m_t =\beta_0 + \beta_1t + \beta_2t^2 + \beta_3t^3\). Entonces podríamos ajustar \(m_t\) mediante mínimos cuadrados ordinarios.

Una extensión de la regresión polinomial es dividir primero el tiempo \(t = 1,. . . , n\), en k intervalos, \([t_0 = 1, t_1]\), \([t_1 + 1, t_2],\cdots,\) \([t_{k − 1} + 1, t_k = n]\); los valores \(t_0\), \(t_1\), …, \(t_k\) se llaman nodos. Luego, en cada intervalo, se ajusta una regresión polinomial, normalmente de orden 3, y esto se llama splines cúbicos. Un método relacionado es suavizar splines, que minimiza el compromiso entre el ajuste y el grado de suavidad dado por

\[\sum_{t=1}^{n}[x_t-m_t]^2+\lambda\int(m_t'')^2 dt,\] donde \(m_t\) es un spline cúbico con nodos en cada tiempo t y el grado de suavidad es controlado por \(\lambda>0.\) El parámetro de suavizado en R es controlado por el argumento spar de la función smooth.spline del paquete stats.

plot(soi)
  lines(smooth.spline(time(soi), soi, spar=.5), lwd=2, col=4)
  lines(smooth.spline(time(soi), soi, spar= 1), lty=2, lwd=2, col=2)

La diferencia entre uan y otra gráfica es es grado de suavizamiento.

<> Tarea: Hacer lo mismo para la serie Recruitment.